Exploratory analysis

From the explatory graph we can see that the gender and age distribution is rather skewed in the data. There are almost 50% more female students than males – especially among the students in their twenties. Even though the mean age is 25.51, there are lots of older students too, the maximum age being 55.

Visual inspection indicates that even though there are slight differences in the attitude and learning approch scores (deep, surface, strategic) between males and females it seems that these differences do not have an effect on median test points. From the figure below we can see that the bulk of the data sits right in the middle of the scale so that there are very few datapoints in the ends of the scales.

Regression analysis

In building linear regression models, it is desirable to have high correlations between the prediction covariates and the response variable, but small correlations between the different prediction covariates. Large correlations between prediction variables leads to the problem of collinearity in linear regression, eg. the model might overfit the data.

Based on explorative analysis in Fig 1, three variables were chosen as prediction covariates.

For the first model, the following variables were chosen as predictors.

• attitude (correlation with points, 0.437)
• strategic learning approach (correlation with points, 0.146)
• surface learning approach (correlation with points, -0.144)

Attitude and strategic learning approach do not correlate (0.06), strategic learning and surface learning approach correlate slightly (-0.16) as well as surfarce learning approach and attitude (-0.17).

model_attitude_stra_surf <- lm(
  data = learning2014, 
  formula = points ~ attitude + stra + surf )

model_attitude_stra <- lm(
  data = learning2014, 
  formula = points ~ attitude)

Regression diagnostics

Regression diagnostics plots help to observe, if the model meets the assumptions of linear regression:

• Linearity of residuals
• Independence of residuals
• Normal distribution of residuals
• Equal variance of residuals

Residuals vs Fitted -plot shows if residuals have non-linear patterns. If the residuals have constant variance of errors, the plot scatters randomly without any distict patterns. We don’t see patterns here, so we can conclude that the variance is constant and the assumption of equal variance of residuals is met.

Normal Q-Q -plot shows if residuals are normally distributed. If the residuals deviate from the straight line it is an indication of non-normality. There are few outlying values in both ends of the line but majority of of the points line nicely and therefore we can conclude that the the residuals of the model are normally distributed.

Residual vs. leverage -plot helps to find influential cases. Not all outliers have leverage in linear regression analysis. Values to observe are at the upper right corner or at the lower right corner – they are the ones which can have an effect on regression results. We don’t see outlying values in these critical spots, so we can conclude that there is no influencial individual cases in the data.

Extra - tweaking the model

Even though we just concluded that there’s no influencial individual cases in the data, we can see from the plots that there are few outlying cases in non-critical spots of the Residual vs. leverage plot. What happens if we exclude them from the data? To do it, we have to single out the cases from the data using Cook’s distance.

Just for the sake of curiosity, let’s investige how excluding three outlier cases (36, 56 & 71) has an effect on our regression model.

learning2014$id <- 1:dim(learning2014)[1]
model_attitude_stra_excluded <- lm(
  data = learning2014 %>% filter(!id %in% c(35, 56, 71)), 
  formula = points ~ attitude)

For the model created using data without cases 35, 56, and 71, the F Test summary: (R²=0.2555, F(1,161)=55.26, p<1e-05). Adjusted R² is 0.25.

The fit of the model seems to be marginally better when three outlying cases are excluded. Let’s see from the data what kind of cases these are so that we could decide if it is prudent to consider them outliers.

Looking the data, two unfortunate students have scored very low in the test, even though their motivation is strong. It might be that some external cause has affected their test performance as other low-scoring (points <=10) students somewhat differ from these cases. There is also one student whose attitude score is very low but who has achieved high points in tests.

Modeling

## 
## Call:
## glm(formula = high_use ~ goout + sex + absnt + studytime, family = "binomial", 
##     data = alc)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0011  -0.7391  -0.4707   0.7279   2.6452  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -2.34971    0.69745  -3.369 0.000754 ***
## gooutlow         0.39022    0.70228   0.556 0.578449    
## gooutmod         0.63488    0.69198   0.917 0.358895    
## goouthigh        2.16469    0.69593   3.110 0.001868 ** 
## gooutvhigh       2.38300    0.71257   3.344 0.000825 ***
## sexM             0.76948    0.27549   2.793 0.005220 ** 
## absnt2           0.41615    0.41509   1.003 0.316074    
## absnt3–4         0.03195    0.38525   0.083 0.933897    
## absnt5–7        -0.20080    0.42637  -0.471 0.637679    
## absnt>8          1.05424    0.36589   2.881 0.003961 ** 
## studytime2–5 h  -0.44563    0.30751  -1.449 0.147302    
## studytime5–10 h -0.91718    0.48515  -1.891 0.058689 .  
## studytime>10 h  -1.01073    0.62811  -1.609 0.107583    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 465.68  on 381  degrees of freedom
## Residual deviance: 374.51  on 369  degrees of freedom
## AIC: 400.51
## 
## Number of Fisher Scoring iterations: 4

alc_glm_second <- glm(
  formula = high_use ~ goout + sex + absnt,
  family = "binomial", 
  data = alc)
summary(alc_glm_second)

## 
## Call:
## glm(formula = high_use ~ goout + sex + absnt, family = "binomial", 
##     data = alc)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9016  -0.7054  -0.4885   0.6704   2.4222  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.72000    0.67394  -4.036 5.44e-05 ***
## gooutlow     0.20910    0.69706   0.300 0.764194    
## gooutmod     0.42097    0.68290   0.616 0.537598    
## goouthigh    1.99006    0.68527   2.904 0.003684 ** 
## gooutvhigh   2.24070    0.70637   3.172 0.001513 ** 
## sexM         0.94468    0.26159   3.611 0.000305 ***
## absnt2       0.44525    0.41055   1.085 0.278133    
## absnt3–4     0.09012    0.37728   0.239 0.811208    
## absnt5–7    -0.15888    0.41776  -0.380 0.703710    
## absnt>8      1.16366    0.36153   3.219 0.001288 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 465.68  on 381  degrees of freedom
## Residual deviance: 379.98  on 372  degrees of freedom
## AIC: 399.98
## 
## Number of Fisher Scoring iterations: 4

Performance comparison

loss_func <- function(class, prob) {
  n_wrong <- abs(class - prob) > 0.5
  mean(n_wrong)
}

training_error <- loss_func(class = alc$hu, prob = alc$predict_high_use)
training_error_simple <- loss_func(class = alc_simple$hu, prob = alc_simple$predict_high_use)

The model has 0.212 training error. Error for initial model is 0.204 so it seems that keeping studytime predictor in the model would have been useful.

library(boot)
cv <- cv.glm(data = alc, cost = loss_func, glmfit = alc_glm_second, K = 10)
cv_simple <- cv.glm(data = alc_simple, cost = loss_func, glmfit = alc_glm_first, K = 10)

Based on the 10-fold cross-validation, the prediction error for the model is 0.22 so it seems that this model performs a bit better than the model introduced in the DataCamp which had 0.26 error. Surprisingly, prediction error for the initial model is 0.22, so it is the best model of all three. To conlude, this somewhat unexpected result idicates that it might be useful to keep also marginally signitficant predictors in the model, as was the case with studytime predictor.

References

Exploring the data

Let’s have a closer look of variables which have higehst correlation with criminality.

Looking the data my intuition is that low-status (likely working-class) housing near highways and industrial areas are more prone to have a higher crime rate.

# Create a categorical variable of the crime rate in the
# Boston dataset (from the scaled crime rate). Use the quantiles as the break
# points in the categorical variable. Drop the old crime rate variable from the
# dataset. Divide the dataset to train and test sets, so that 80% of the data
# belongs to the train set. (0-2 points)

boston_scaled %<>% 
    mutate(
      crime_cat = quantcut(
        per_capita_crime, 
        labels = c("very low","low","high","very high")
        )
      ) %>% select(-per_capita_crime)


# Let's use  Caret to create training dataset and make sure that
# the overall class distribution is preserved
# eg. the random sampling occurs within each class

rows_to_train <- caret::createDataPartition(boston_scaled$crime_cat, p = .8, 
                                  list = FALSE, 
                                  times = 1)

boston_training <- boston_scaled[rows_to_train,]
boston_test <- boston_scaled[-rows_to_train,]

Linear discriminant analysis and prediction

The model predicts the correct class with 73 % accuracy. Location near to highway entraces seems to contribute most to the criminality.

K-means clustering

Plot represents the variance within the clusters. It decreases as k increases, but one can notice a bend (or “elbow”) at ~2. This bend indicates that additional clusters beyond this point have only little value. So, let’s fit k-means clustering with k=2.

Looking the plots, our two clusters correspond pretty nicely to our intuition of crime-prone neighbourhoods residing near high-pollution industrial areas and highways. More affluent cluster 1 (colour black) has less criminality. Cluster 2 has more pollution, is located near to highways, has factories nearby and has more criminality.

References

Exploring human development data with principal component analysis

The dataset we’ll use is from the United Nations Development Programme. For the following analysis, the dataset was preprocessed using an external script.

The dataset has the following variables:

Most of the variable distributions are more-or-less skewed, so it is very likely that PCA will not work without normalizing the data. Only the edu_exp (expected years of education) variable looks normally distributed.

Unexpectedly, maternal mortality and life expectency are highly correlated (-0.86) as well as high birth rate and maternal mortality (0.76). Variables regarding education have also strong correlation with maternal mortality, so that higher values in education variables, lead to lower value in maternal mortality rate. Secondary education of women have an association to lower birth rate (-0.68), as well as higher life expectancy (-0.72). Quite surprisingly, women’s representation in the parliament doesn’t seem to have an association with any of other variables or the association is very low. As one might expect, maternal mortality and high birth rate are associated with low GNI.

pca_summary_rows <- c(
  "Standard deviation",
  "Eigenvalue",
  "Proportion of Variance",
  "Cumulative Proportion")

pca_human <- prcomp(human)
pca_summary <- pca(pca_human)
rownames(pca_summary) <- pca_summary_rows

# Use built-in scaling in prcomp() -function
pca_human_final <- prcomp(human, scale = TRUE, center = TRUE)
pca_summary_final <- pca(pca_human_final)
rownames(pca_summary_final) <- pca_summary_rows

Exploring tea survey data with multiple correspondence analysis

The following questions which were used to gather the tea dataset can be found from the book Exploratory Multivariate Analysis by Example Using R, Second Edition, pages 131–132.

1. What kind of tea do you drink the most (black tea, green tea, flavored tea)
2. How do you take your tea (nothing added, with lemon, with milk, other)
3. What kind of tea you buy ? (tea bags, loose tea, both)
4. Do you add sugar to your tea (yes, no)?
5. Where do you buy your tea (supermarket, specialist shop, both)
6. What kind of tea you buy (cheapest supermarket brand, well-known brand, upscale, it varies)
7. How ofter do you drink tea (more than twice a day, once a day, 3 to 6 times a week, once or twice per week)?

Location in wich the tea is drunk (yes/no):

8. Do you drink tea at home?
9. Do you drink tea at work?
10. Do you drink tea in tearooms or coffee shops?
11. Do you drink tea at frinds houses?
12. Do you drink tea in restaurants?
13. Do you drink tea in bars?

Time of the day (yes/no):

14. Do you drink tea at breakfast?
15. Do you drink tea in the afternoon?
16. Do you drink tea in the evening?
17. Do you drink tea after lunch?
18. Do you drink tea after dinner?
19. Do you drink tea throughout the day?

Questions about the image of the product (yes/no):

20 Do you consider tea to be exotic? 21. Do you associate tea with spirituality? 22. Is tea good for your health? 23. Is tea a diuretic? 24. Do you associate tea with friendliness? 25. Does tea stop the body from absorbing iron? 26. Is tea feminine? 28. Is tea refined? 29. Will tea help you to lose weight? 30. Is tea a stimulant? 31. Is tea a relaxant? 32. Does tea have any effect on your overall health?

Additionally, there were four background questions regarding:

• sex
• professional category (farmer, manual labourer, professional, senior management, employee, other profession, student, unemployed)
• age
• regular participation in sports (yes or no).

Skim summary statistics
n obs: 300
n variables: 36

Let’s use Goodman-Kruskal tau to explore associations between variables in the dataset. Goodman-Kruskal tau can be used to study association when the measurement level of the variable is at least categorical.

Going trough the Goodman-Kruskal tau-matrix, we can observe few interesting associations between variables. Overall, associations are not very strong, but the following associations emerge when taking into account associations which are stronger than the mean tau-value (0.08).

• gender and considering tea feminine or not feminine
• age and considering tea feminine or not feminine
• age and a belief that tea has an effect on overall health
• a belief that tea has an effect on overall health and and considering tea healthy
• type of the tea and age
• type of the tea and habit of using sugar
• drinking tea at tearoom and buying tea from specialist shops
• frequency of drinking tea and drinking tea at breakfast
• price, type of tea and the place of purchase

Even though we’ll eventually use all variables in the MCA, let’s observe few of these associations bit more closely, just for the sake of curiosity.

• There are more female than male responders to the survey. Over the half of the responders are under 35.
• It seems that the bulk of the tea is being bought from the chain stores and the preferred style is Earl Gray.
• Interestingly, considering all answers tea is not considered feminine, but females tend to associate tea with femininity more often than men.

tea_mca <- MCA(tea,
    quanti.sup = 19, # age as quantitative supplement
    quali.sup = c(20:36),
    graph = FALSE) # other background variables as categorical 
                                # supplements

ggord(tea_mca,
      tea$Tea,
      arrow = 0.2,
      txt = 3,
      vec_ext = 1.05,
      size = 2,
      repel = TRUE) + theme_minimal()